En el flujo con una cavidad rectangular, el fluido se encuentra confinado en una región rectangular de dimensiones $L\times W$ donde el movimiento se debe al desplazamiento transversal de una de las caras. En este sentido, la velocidad de la cara superior es de $U_0$ y el de las demás interfaces es de $0$. De acuerdo con esto, se tienen las siguientes consideraciones:
De aquí, las ecuaciones de continuidad y de Navier-Stokes se reducen a:
\begin{align} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho \vec{u})=0 \quad &\to\quad \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 \\ \rho\left(\frac{\partial \vec{u}}{\partial t}+\vec{u}\cdot\nabla\vec{u}\right)=\rho\vec{f}-\nabla p+\mu\nabla^2\vec{u}\quad &\to\quad \begin{cases} \rho \left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) \\ \rho \left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{\partial p}{\partial v}+\mu\left(\frac{\partial^2v}{\partial x^2}+\frac{\partial^2v}{\partial y^2}\right) \end{cases} \end{align}Se puede ver entonces que la ecuación resultante es un caso particular del problema de convección-difusión estable:
$$\frac{\partial f}{\partial t} + \vec{u}\cdot\nabla f=k\nabla^2f+\frac{s}{\rho c_p}$$Donde $\frac{\partial f}{\partial t}=0$, $k = \frac{\mu}{\rho}=\nu$, y $\frac{s}{c_p}=-\frac{\partial p}{\partial *}$ para $f=u, *=x$ y $f=v, *=y$, y las condiciones de frontera corresponden a valores conocidos (tipo Dirichlet).
Ahora, para aplicar el método de elementos espectrales, es necesario realizar los siguientes pasos:
Cada uno de los pasos se tratará en detalle a continuación
In [1]:
# Prepara el ambiente de trabajo con las librerías a utilizar
import numpy as np # Ofrece funciones para operaciones básicas con arreglos
import scipy.linalg as lin # Ofrece funciones para operaciones de álgebra lineal
import matplotlib.pyplot as plt # Permite incluir un ambiente de visualización
from mpl_toolkits.mplot3d.axes3d import Axes3D
import timeit # Permite la medición de tiempos para los distintos algoritmos
from sem2D import * # Agrupa las funciones externas en un archivo
# from __future__ import division # Corrige posibles problemas que surgen a la hora de dividir enteros
#%pylab inline
In [2]:
## Valores de prueba
L = 2.0
W = 2.0
mu = 1.
U0 = 1.
In [3]:
Ne = np.array([12,8]) # Indica la cantidad de elementos a lo largo de los ejes x y y
Np = np.array([3,2]) # Indica la cantidad de divisiones sobre cada elemento en xi y en eta
# Genera los puntos correspondientes a los vertices de cada elemento
(px, hx) = np.linspace(-L/2,L/2,Ne[0]+1,retstep=True)
(py, hy) = np.linspace(-W/2,W/2,Ne[1]+1,retstep=True)
Conociendo el número de elementos y su distribución, es posible generar los nodos correspondientes a los vértices de cada elemento. Sin embargo, como no se trabaja con elementos lineales, es necesario generar los nodos de interpolación para cada elemento. Para ello, es necesario definir la familia de polinomios sobre la cual van a estar basados estos nodos. En este caso se utilizan los nodos de interpolación de Lobatto.
In [4]:
(x, y, ne, npe) = genNodes(Ne,Np,px,py)
# Grafica la malla de puntos
for l in range(ne):
plt.plot([x[l,0],x[l,Np[0]],x[l,-1],x[l,-Np[0]-1],x[l,0]],
[y[l,0],y[l,Np[0]],y[l,-1],y[l,-Np[0]-1],y[l,0]],
'r--')
plt.hold(True)
plt.plot(x,y,'bo',ms=4)
plt.xlabel('coordenada x')
plt.ylabel('coordenada y')
plt.show()
In [5]:
h = min(hx,hy)
(coords, C, gfl, ng) = genConn(Ne,Np,ne,npe,x,y,L,W,U0,h,'quad')
Para la generación de la matriz del problema es necesario definir las matrices de difusión, masa y advección por cada elemento:
$$A_{ij}^{(\ell)}=\iint_{E_\ell} \nabla\psi_i^{(\ell)}\cdot\nabla\psi_j^{(\ell)}dxdy \quad B_{ij}^{(\ell)}=\iint_{E_\ell} \psi_i^{(\ell)}\psi_j^{(\ell)}dxdy \quad C_{ij}^{(\ell)}=\iint_{E_\ell} \psi_i^{(\ell)}\vec{u}\cdot\nabla\psi_j^{(\ell)}dxdy$$Estas integrales se obtienen utilizando la cuadratura de Lobatto sobre cada una de las dimensiones, ya que las funciones de interpolación son bilineales y se escriben de la forma:
$$ \psi_i(\xi,\eta) = \mathcal{L}_i(\xi)\mathcal{W}_i(\eta) $$donde $\mathcal{L}_i(\xi)$ y $\mathcal{W}_i(\eta)$ son las funciones de interpolación (unidimensionales) correspondientes. Por simplicidad, se trabaja sobre el elemento estandar con coordenadas $(-1,-1),(1,-1),(-1,1) \text{ y }(1,1)$ y se alude a un cambio de coordenadas.
Comenzamos con la matriz que define el sistema de ecuaciones, o la matriz de difusión. Dado que esta depende de los gradientes de las funciones de interpolación, es necesario realizar todos los pasos de construcción de los términos como sigue.
En este caso, debido al orden variable del método aplicado, se generan las funciones de interpolación a partir de los polinomios de interpolación de Lagrange. Para simplificar la notación se considera la biyección: $\{0,1,\dots,Np\} \to \{0,\dots,m\}\times\{0,\dots,n\}$ definida por $i \mapsto (j,k)$
$$ \psi_i(\xi,\eta) = \psi_{(j,k)}(\xi,\eta) = \mathcal{L}_j(\xi)\mathcal{W}_k(\eta) = \left(\prod_{l\neq j}^m\frac{\xi-\xi_l}{\xi_j-\xi_l}\right) \left(\prod_{l\neq k}^n\frac{\eta-\eta_l}{\eta_k-\eta_l}\right)$$De acuerdo con esta descripción, se tiene entonces que las derivadas con respecto a $\xi$ y $\eta$ son entonces:
\begin{align} \frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &= \psi_{(j,k)}(\xi,\eta) \sum_{l\neq j}^m \frac{1}{\xi-\xi_l} = \left[\sum_{l\neq j}^m \left(\frac{1}{\xi_j-\xi_l}\prod_{p\neq j,l}^m \frac{\xi-\xi_p}{\xi_j-\xi_p}\right)\right]\left(\prod_{l\neq k}^n\frac{\eta-\eta_l}{\eta_k-\eta_l}\right) \\ \frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} &= \psi_{(j,k)}(\xi,\eta) \sum_{l\neq k}^m \frac{1}{\eta-\eta_l} = \left(\prod_{l\neq j}^n\frac{\xi-\xi_l}{\xi_j-\xi_l}\right) \left[\sum_{l\neq k}^m \left(\frac{1}{\eta_k-\eta_l}\prod_{p\neq k,l}^m \frac{\eta-\eta_p}{\eta_k-\eta_p}\right)\right] \end{align}De este modo, la función de transformación de coordenadas corresponde a:
$$ \vec{x} = \sum_{i=0}^{N_p-1} \vec{x}_i^E \psi_i(\xi,\eta) = \sum_{i=0}^{N_p-1} \vec{x}_i^E \psi_{(j,k)}(\xi,\eta) = \sum_{i=0}^{N_p-1} \vec{x}_i^E \mathcal{L}_j(\xi)\mathcal{W}_k(\eta) $$y por lo tanto sus derivadas son:
\begin{align} \frac{\partial x}{\partial \xi} &= \sum_{i=0}^{N_p-1}x_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &\quad \frac{\partial x}{\partial \eta} &= \sum_{i=0}^{N_p-1}x_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} \\ \frac{\partial y}{\partial \xi} &= \sum_{i=0}^{N_p-1}y_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\xi} &\quad \frac{\partial y}{\partial \eta} &= \sum_{i=0}^{N_p-1}y_i^E\frac{\partial\psi_{(j,k)}(\xi,\eta)}{\partial\eta} \end{align}Ya con estos valores es posible soluciones los gradientes de las funciones solucionando los sistemas:
$$ \begin{pmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{pmatrix} \nabla\psi_i = {\frac{\partial\psi_i}{\partial\xi} \choose \frac{\partial\psi_i}{\partial\xi}} $$De manera que el coeficiente de cambio de coordenadas corresponde a:
$$ Det[J] = \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{vmatrix} = \frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta}-\frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \xi} $$Para el otro lado de la ecuación se implementa la función de entrada. Para esto se recorre sobre todos los elementos y se van calculando los valores para la función de manera correspondiente. Estos valores se acomodan al sistema global aprovechando la matriz de conectividad.
Así mismo, es necesario modificar la matriz global de manera que se incorporen las condiciones de frontera de Dirichlet. Si $l$ es un nodo de frontera, esto se hace a través de los pasos:
In [6]:
(Mat_dif, vec_b) = matDiff(Ne,Np,ne,npe,ng,x,y,C,gfl,mu)
plt.spy(Mat_dif)
plt.show()
Teniendo ya construidos los componentes del sistema a resolver, basta aplicar un esquema de solución razonable para el tipo de matriz obtenida. En este caso se utiliza el solucionador por defecto del módulo linalg. Dado que la matriz construida inicialmente es singular, se omiten los valores de la primera fila y columna - que corresponden a un nodo de frontera y ya se conoce la velocidad en este punto. Esto permite que la matriz obtenida sea invertible.
In [7]:
start = timeit.default_timer()
sols = lin.solve(Mat_dif[1:,1:],vec_b[1:])
print timeit.default_timer()-start
# Recupera las velocidades
solx = np.concatenate([[gfl[0,1]], sols[0:ng-1]])
soly = sols[ng-1:2*ng-1]
# Calcula la magnitud de la velocidad en cada nodo
solm = np.abs(solx**2+soly**2)
# Recupera las presiones
solp = sols[2*ng-1:]
In [8]:
# Grafica las fronteras de los elementos
for l in range(ne):
plt.plot([x[l,0],x[l,Np[0]],x[l,-1],x[l,-Np[0]-1],x[l,0]],
[y[l,0],y[l,Np[0]],y[l,-1],y[l,-Np[0]-1],y[l,0]],
'r--')
plt.hold(True)
# Grafica la magnitud del campo de velocidades
for l in range(ne):
Px = (coords[C[l,:],0]).reshape(Np[1]+1,npe/(Np[1]+1))
Py = (coords[C[l,:],1]).reshape(Np[1]+1,npe/(Np[1]+1))
SM = (solm[C[l,:]]).reshape(Np[1]+1,npe/(Np[1]+1))
plt.pcolor(Px,Py,SM,vmin=min(solm), vmax=max(solm))
# Grafica el campo de velocidades
plt.quiver(coords[:,0],coords[:,1],solx,soly,color='w')
plt.title('Campo de Velocidades para el Flujo en la Cavidad')
plt.xlabel('Posicion (coord x)')
plt.ylabel('Posicion (coord y)')
plt.xlim([-1., 1.])
plt.ylim([-1., 1.])
plt.show()
In [9]:
# Grafica el campo de presiones
xpres = np.zeros(ne)
ypres = np.zeros(ne)
for l in range(ne):
xpres[l] = (x[l,0]+x[l,Np[1]]+x[l,-1]+x[l,-Np[1]-1])/4
ypres[l] = (y[l,0]+y[l,Np[1]]+y[l,-1]+y[l,-Np[1]-1])/4
xpres = xpres.reshape(Ne[1],Ne[0])
ypres = ypres.reshape(Ne[1],Ne[0])
zpres = solp.reshape(Ne[1],Ne[0])
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(xpres,ypres,zpres)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Presion')
#plt.pcolor(xpres,ypres,zpres)
plt.show()
In [ ]: